path(path,'FonctionsFormes');
path(path,'Dirichlet');
path(path,'Erdc');
path(path,'FluxFree');
path(path,'PostProcess');
path(path,'ZZ');
path(path,'Globale');
path(path,'Remaillage');
path(path,'Solve');
path(path,'ErreurLocale');


% Données géométriques
%%%%%%%%%%%%%%%%%%%%%%

% Poutre en acier, longueur L, section circulaire de 40cm de diamèetre

% Longueur (en m)
donnees.poutre.L = 10 ;

% Propriétés de la poutre
donnees.poutre.E = 210 * 1000000000 ;
donnees.poutre.I = 0.0012 ;

%% Boucle sur le nombre d'éléments

nmin = 3 ;
nmax = 50 ;

efficaciteExplicite = zeros(nmax-nmin+1,1) ;
efficaciteDirichlet = zeros(nmax-nmin+1,1) ;
efficaciteZZ = zeros(nmax-nmin+1,1) ;

tablErrExplicite = zeros(nmax-nmin+1,1) ;
tablErrDirichlet = zeros(nmax-nmin+1,1) ;
tablErrZZ = zeros(nmax-nmin+1,1) ;
tablErrRdC = zeros(nmax-nmin+1,1) ;

tablNbElem = zeros(nmax-nmin+1,1) ;


for n=nmin:nmax
    
    % Données problème
    %%%%%%%%%%%%%%%%%%%
    
    % Type de chargement
    donnees.chargement='ponctuel'; % 'ponctuel' ou 'reparti' ou 'ponctuelETreparti'
    
    % Force (en N)
    donnees.F = 500;
    
    % Force linéique de la forme ax² + bx + c
    donnees.forceLineiqueA = 5 ;
    donnees.forceLineiqueB = 20 ;
    donnees.forceLineiqueC = 10 ;
    
    % Interpolation avec fonctions de forme: (nombre de points d'interpolation
    % par éléments)
    donnees.nInt = 100;
    
    % Nombre d'éléments initial
    donnees.nElemIni = n ;
    tablNbElem(n-nmin+1) = donnees.nElemIni ;
    donnees.nombreElements = donnees.nElemIni ;
    
    % Nombre de subdivision d'un élément pour la resolution par patchs
    donnees.subdivisionPatch = 10 ;
    
    % Discretisation de la solution analytique
    donnees.discretisationAnal = 100;
    
    % Maillage du domaine
    donnees.maillage = Mailleur_uniforme(donnees) ;
    maillage = donnees.maillage ;
    
    % Création de la matrice de rigidité de de la matrice de contrainte (prise
    % en compte des CL par lagrangien
    K = AssemblageRigidite(maillage,donnees.poutre) ;
    
    C = [1,zeros(1,2*length(maillage)-1);
        0,1,zeros(1,2*length(maillage)-2)];
    
    Klag = [K ,C';C,[0,0;0,0]];
    
    force = [zeros(2*(length(maillage)-1),1);donnees.F;0] ;
    force = force + second_membre_lineique(donnees,maillage) ;
    
    forceLag = [force;0;0];
    
    U = Klag\forceLag ;
    donnees.res.Uh = U(1:length(U)-2) ;
    %Usubs = Ksubs\forceSubs ;
    
    % plot(donnees.res.Uh);
    
    %derivSecUhL(donnees)
    %fleche = U(length(U)-3)
    
    
    
    [donnees.UhInt,donnees.maillageInt] = Fonction_interpolee(donnees);
    
%     Uanal = resultatAnalytique(donnees);
%     flecheAnal = Uanal(length(Uanal))
    
    [errRdC errRdCelem] = ERdC(donnees) ;
    tablErrRdC(n-nmin+1) = errRdC ;
    
    erreurResiduExplicite = EvalErrGlob(donnees) ;
    tablErrExplicite(n-nmin+1) = erreurResiduExplicite ;
    efficaciteExplicite(n-nmin+1) = (erreurResiduExplicite)/errRdC ;
    
    [erreurPatchDirichlet, erreurPatchDirichletElem] = Err_globale_patch_Dirichlet(donnees) ;
    tablErrDirichlet(n-nmin+1) = erreurPatchDirichlet ;
    efficaciteDirichlet(n-nmin+1) = (erreurPatchDirichlet)/errRdC ;    
    
    [errZZ errZZelem Mzz] = EZZ(donnees) ;
    tablErrZZ(n-nmin+1) = errZZ ;
    efficaciteZZ(n-nmin+1) = (errZZ)/errRdC ; 
    
%     Moments_interpoles(donnees) ;
%     hold on ;
%     plot(donnees.maillage,Mzz) ;
    
end

figure;
plot(tablNbElem,efficaciteExplicite);
title('efficacité erreur explicite');

figure;
plot(tablNbElem,efficaciteDirichlet);
title('efficacité erreur patch Dirichlet');

figure;
plot(tablNbElem,efficaciteZZ);
title('efficacité erreur ZZ');

figure;
plot(tablNbElem,efficaciteExplicite);
title('efficacité erreur explicite');

figure;
plot(tablNbElem,tablErrRdC);
title('Evolution de l erreur exacte en fonction de la taille des éléments');

figure;
plot(tablNbElem,tablErrExplicite);
title('Evolution de l erreur explicite en fonction de la taille des éléments');

figure;
plot(tablNbElem,tablErrDirichlet);
title('Evolution de l erreur par patch de Dirichlet en fonction de la taille des éléments');

figure;
plot(tablNbElem,tablErrZZ);
title('Evolution de l erreur ZZ en fonction de la taille des éléments');



%plot(donnees.maillageInt,donnees.UhInt)
%plot(1:donnees.nombreElements, errRdCelem)
%hold on ;
%plot(0: donnees.poutre.L/(donnees.discretisationAnal) :donnees.poutre.L ,Uanal,'-r')


